calcBalancePrecision <- function(df, te = 1, sd = 1, coefs = c(1,2,3)){
    
  ## store Hansen Bowers 08 d^2 p-value, Atkinson 03 hz balance and var(delta-hat) precision, and RMSE
  
  sig2 <- sd^2

  cov.names <- names(df)[grep("cont", names(df))]
  x <- as.matrix(cbind(1, df[, cov.names]))
  units <- nrow(x)
  
  t <- as.numeric(as.factor(df$Tr))
  y <- cbind(x, t)%*%c(coefs, te)+rnorm(units, 0, sd)
  lm.out <- lm(y ~ cont1 + cont2 + t, data = data.frame(cbind(x, t)))
  
  ## vdh
  a <- rep(NA, units)
  a[df$Tr=="Treatment 1"] <- 1
  a[df$Tr=="Treatment 2"] <- -1
  vdh <- sig2*solve(a%*%a - t(a)%*% x%*% solve(t(x)%*%x)%*% t(x)%*% a)
  
  ## hz 
  H <- matrix(0, units, 2)
  for(i in 1:units){
    if(df$Tr[i]=="Treatment 1"){
      H[i,1] <- 1
    }else{
      H[i,2]<-1
    }
  }
  hz <- t(H)%*%x[, cov.names]   
  hz <- det(hz)
  
  ## d2.p
  ffff <- as.formula(paste("as.numeric(as.factor(Tr)) ~ ", cov.names, sep = ""))
  xb <- xBalance(ffff, data = df, report=c("chisquare.test"))
  d2.p <- xb$overall$p.value
    
  ## rmse
  ## difference in means:
  mean.y2.obs <- mean(y[df$Tr == "Treatment 2"])
  mean.y1.obs <- mean(y[df$Tr == "Treatment 1"])
  rmse <- sqrt(mean(c((y[df$Tr == "Treatment 2"]-mean.y2.obs), (y[df$Tr == "Treatment 1"]-mean.y1.obs))^2))
  ## lm:
  #  rmse <- sqrt(mean(c((y[df$Tr == "Treatment 2"] - lm.out$fitted.values[df$Tr == "Treatment 2"]), (y[df$Tr == "Treatment 1"] - lm.out$fitted.values[df$Tr == "Treatment 1"]))))
  
  return(data.frame(d2.p = d2.p, hz = hz, vdh = vdh, rmse = rmse))
}
